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We theoretically study the scattering of a plane wave of a 
ballistic electron on a circular n-p junction in single and 
bilayer graphene. We compare the exact wave function 
inside the junction to that obtained from a semiclassi- 
cal formula developed in catastrophe optics. In the semi- 
classical picture short-wavelength electrons are treated 
as rays of particles that can get reflected and refracted at 
the n-p junction according to Snell's law with negative 
refraction index. 



We show that for short wavelength and close to caus- 
tics this semiclassical approximation gives good agree- 
ment with the exact results in the case of single-layer 
graphene. We also verify the universal scaling laws that 
govern the shrinking rate and intensity divergence of 
caustics in the semiclassical limit. It is straightforward 
to generalize our semiclassical method to more com- 
plex geometries, offering a way to efficiently design and 
model graphene-based electron-optical systems. 
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1 Introduction The curves of caustics, which are 
envelopes of a family of rays at which the density of rays is 
singular, can be described using geometrical optics. How- 
ever, the wave intensity in the vicinity of caustics cannot be 
predicted by the simplest theory where rays are endowed 
with amplitude and phase, and allowed to interfere where 
they cross: as shown by Berry and Upstill lUl, such cal- 
culations fail exactly on caustics by predicting intensity 
divergences. They have shown that applying catastrophe 
theory ||2l in optical systems can offer a good solution: 
catastrophe optics predicts finite intensity on the caustics 
and gives quantitative results in the short-wavelength (or 
'semiclassical') limit. This tool enables the classification 
of caustics according to their structural stability. Each class 
has its own characteristic diffraction pattern, and the cor- 
responding wave function can be well approximated by an 
integral representation in terms of a polynomial describ- 
ing the class 1 1 1. Near the caustics, the so-called universal 
scaling laws govern the /c ^ oo asymptotic behavior of the 
wave function, where k is the wave number. In particular, 
as k increases the intensity also increases, and the diffrac- 



tion fringe dimensions decrease proportionally to certain 
power functions of k with universal exponents. 

Many theoretical ||3 and experimental |4| works have 
already been published on materials with negative re- 
fractive index. Recently, it has been demonstrated that in 
graphene-nanostructures, the optics of electron flow can 
be described by reflections and refractions with negative 
refractive index |5|, moreover, it has also been shown that 
in such systems, caustics can arise |[5l[6l|7l. 

With the example of a circular n-p junction (NPJ) (61 
[71, we investigate in this paper whether catastrophe optics 
could be applied in graphene in the short wavelength limit. 
Our main finding is that in single-layer graphene (SLG), 
near the caustics the wave function can be well approxi- 
mated using an integral formula developed in catastrophe 
optics. This offers an efficient way to design and model 
graphene-based electron-optical devices. 

2 Catastrophe theory and integral representation 
of the wave function Let us consider the scattering of 
ballistic electrons in graphene in the setup shown in Fig.[Tl 
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Figure 1 An incident 
plane wave of a ballistic 
electron is scattered by 
a circular potential bar- 
rier V(r). The incoming 
electron on the n side is 
in the conduction band, 
while on the p side in the 
valence band. 



We use a circular NPJ, formed by a simple gate potential 
V{r) = VoO{R — r), where O denotes the Heaviside func- 
tion, R is the radius of the junction and the energy E of the 
incident particles fulfills < E < Vq. The assumption 
of such a sharp potential is valid provided that Xf ^ d, 
where Xf is the Fermi wavelength and d is the character- 
istic length scale in which the scattering potential varies. 
To prevent intervalley scattering, we assume that d ^ a, 
where a is the lattice constant of graphene (SUS) . The ex- 
act calculations of the wave function inside and around the 
junction are presented in Ref. |6J for SLG and in Ref. [7] 
for bilayer graphene (BLG). In the same papers it has been 
shown that caustics clearly arise in the wave function pat- 
tern. This suggests that we consider the propagation of 
electrons as that of rays of light, and caustics as the en- 
velopes of these rays. These can be refracted and reflected 
according to Snell's law: 



sma 
sin/3 



/Cij 



n, 



(1) 



where a and P are the angles of incidence and refraction, 
respectively, and /cin and k are the wave numbers inside 
and outside the junction. Just like in certain photonic crys- 
tals ||9l, in graphene NPJs, the refractive index n is nega- 
tive as shown in Ref. . The curves of the caustics can be 
calculated using differential geometry, and the results are 
given by Eq. (9) in Ref. (6]. Note that the caustics' curves, 
which are the same in SLG and in BLG, depend only on 
n. For any given n, we can label the curves of the caus- 
tics according to the number of internal reflections of the 
rays that create them. The first order caustic (p = 1) cor- 
responds to the envelope of the rays that were not reflected 
inside the junction, and the second order caustic (p = 2) is 
formed by the rays after one reflection. Thus p — 1 gives 
the number of the internal reflections. 

We now consider the problem in terms of catastrophe 
optics. As catastrophe theory O is a special case of the 
more general singularity theory, catastrophe optics deals 
with singularities of ray families called caustics (for details 
see, e.g., Ref. |1|). 

Let us denote the optical length by (/), which is the path 
integral weighted by the position-dependent refraction in- 
dex n: 0(a, 9^) = /i + n/2, for the notation see Fig. [J] To 
find the extrema of the action, we have the following con- 
dition for any given 9^: d(j){a^ 9^) /9a = 0. This results in 
a set of a corresponding to the rays passing through the 




Figure 2 An incoming 
ray of electrons partly 
reflects, and partly en- 
ters the NPJ. The upper 
branch of the p = 1 
fold caustic is shown 
in the 2nd quarter of 
the NPJ. The dashed 
line is normal to the 
junction. d\ denotes the 
observation point. 



point d\. Basically, this is the well-known Fermat's prin- 
ciple. Varying 9^, the extremum points a (9^) will change. 
When two or more extrema coalesce, (j) will be stationary 
to higher than first order. It means 9^ is on a caustic, if and 
only if 9^0(a, 9^)/9a^ = 0. For example, on a fold type 
caustic, two rays touch, while in a cusp type, three. 

For short wavelength, the stationary phase approxima- 
tion can be used, and one can simply sum up the contribu- 
tions of the rays passing through 9^. The problem is that 
this sum diverges exactly on caustics where, however, the 
integral representation (IR) of the wave function behaves 
well, and gives a smooth and accurate solution: 

^l)ir{^) = cVk / da b{a) exp{ik(l){a, 9^)) , (2) 

where b{a) is a weighting function, c is a constant prefactor 
and 7 = arcsin(min(l, |n|)). Although the general form 
of Eq. (|2]) was originally derived from the Helmholtz equa- 
tion im, it can be used in SLG because the SLG eigenstates 
solve the Helmholtz equation. If we neglect the evanescent 
waves, the BLG eigenstates are solutions as well. How- 
ever the evanescent waves cannot be neglected close to the 
boundary of the NPJ, and so we cannot expect that Eq. Q 
would give an accurate result in BLG not even far from the 
boundary. Note that in graphene systems, b{a) is a spinor 
and it can be determined from the boundary conditions fol- 
lowing the derivation presented in Ref. 11]. We determined 
the c prefactor by fitting V^ir(9^) to the exact results. 

The reason for the failure of the stationary phase ap- 
proximation is that the ray contributions to the wave func- 
tion are always considered separately. It is clear from 
Eq. ^ that for the limit k -^ 00, the integrand is a 
rapidly oscillating function of a, and destructive inter- 
ference occurs for every a except for those satisfying 
90(a,9^)/9a = 0. Therefore far from caustics, the sta- 
tionary points are well separated, and the method of sta- 
tionary phase can be used in its simplest form. On the other 
hand, near a caustic, two or more of the stationary points 
lie close together, and as soon as the distance between 
them is comparable to the wavelength of the particles, 
their contributions are not separated anymore by a region 
of destructive interference. In this case, they cannot be 
simply added as in the stationary phase approximation. 
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a) Exact solution: |'0(x, ^)|^ 




b) Semiclassical approximation: |'0ir(x, 2/)|^. 
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Figure 3 Intensity around the focal point in SLG. 
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a) Exact solution: \ip{x,y)\'^ 




b) Semiclassical approximation: \ipir{x,y)\'^. 
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Figure 4 Intensity around the missing focal point in BLG. 



3 Numerical results To check the accuracy of the 
integral approximation Q, we calculate the wave function 
in SLG and in BLG for many different n and p values. In 
each case, we compare the exact results for the probabil- 
ity function |?/;(x, 7/)p inside the NPJ obtained by the ex- 
act calculation (61I3 to the approximated ones. To make a 
quantitative comparison, we introduce the Frobenius norm 

of a matrix A defined as v'Tr(AtA) = JT^^j l^^jp. 

where A^ denotes the adjoint matrix of A. We define the 
error of our approximation as the ratio of the Frobenius 
norm of the difference matrix | ?/^ir (x^ , ?/j ) p — \ijj {xi , ^j ) P 
and that of the matrix | ?/^(xi , ^j ) p , where Xi and i/j are the 
coordinates of the grid points inside the NPJ, |?/^(xi, ?/j)p 
is the exact probability function, and \ijjir{xi^yj)\'^ is the 
electron density obtained from Eq. Q. 

Here we present only the case of p = 1 and n = — 1. 
One can clearly see in Figs. [3] and IH that the approximation 
works well indeed. Here kR = h^R = 10000, the plotted 
range is -0.517 < x < -0.5, -0.0025 < y < 0.0025, 
and the solid line shows the analytic curve O. Note that 
in BLG the focus is missing due to the zero transmission 
probability at perpendicular incidence Q. We calculate 
the error as discussed above, and we find that it is 0.01089 
in SLG, and 0.2038 in BLG. The larger error for BLG is 
due to an apparent deviation from the exact results, namely 
the undulation in Fig.|4al which is missing in Fig.|4bl Note 
that the undulations wavelength is in the order of the parti- 
cles' de Broglie wavelength. 

In Fig. [5] we present results for the fold type caustics. In 
SLG, the intensity of the wave function is calculated along 
a segment parallel to the x-axis passing through a point 
on the fold caustic, at which the tangent of the caustic is 
parallel to the y-Sixis (see inset). According to the theory 
outlined in Ref. d, the wave function along this segment 



should be proportional to the Airy function. Here we inves- 
tigate whether the intensity follows the squared Airy func- 
tion. One can see that the two curves, one fitted on the exact 
and the other on the approximated intensities, are basically 
the same, they cannot be distinguished. 
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Figure 5 Intensity in SLG as the function of x close to the 
fold caustics by kR = h^R = 2500, n = —1, andp = 1. 
Green (bright) points indicate the exact results, and blue 
(dark) ones the approximation. One red (large) point shows 
the crossing point of the segment and the fold caustic at 
(-l/\/2, 1/a/8). The fitted aAi^{bx + c) type functions 
completely overlap for the exact and approximated results 
(solid line). (Here a, 6, and c are the fitting parameters of 
the Airy function.) 
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3.1 Scaling laws Catastrophe theory predicts cer- 
tain scaHng laws that govern the /c ^ oo asymptotic be- 
havior of the wave function, i.e., the shrinking rate of the 
diffraction fringes in the different dimensions and the di- 
vergence of the intensity. Near the caustics, the scaHng 
laws relate the /c-independent standard diffraction catastro- 
phe ^{91) \l\ to the /c-dependent physical wave function 
V^ir(^) given by Eq. Q: 

where r^ are the coordinates set to the caustics in the ith 
dimension, and tz = kR is the dimensionless wave num- 
ber. Here (3 is the singularity index, the measure of the di- 
vergence at the most singular point, where r^ = for all 
i, and ai are related to the fringe spacings in the differ- 
ent control directions. These indices are invariant for each 
classes: /3 = 1/6 and ai = 2/3 for fold catastrophes, while 
/3 = 1/4, cTi = 1/2, and (72 = 3/4 for cusp type caus- 
tics. Note that the shrinking is anisotropic since ai 7^ (72. 
Around the cusp type focal point in the NPJ, the coordi- 
nates (ri, r2) coincide with the coordinates (x, y) with the 
exception that its origin is in the focus. On the fold caustic, 
the only axis of the local coordinate system for the control 
parameter ri is perpendicular to the caustic, and the origin 
is in the crossing point on the caustic. 

We have performed several numerical tests to check the 
scaling laws, here we present two of them. In each cases 
we fitted the theoretical curves on the exactly calculated 
results for the intensity |?/^(x, ^) p. (I) We calculated the in- 
tensity as a function of k in the focus in SLG for p = 1 and 
n = — 1. (Hereafter, the parameter k is in units of l/R.) By 
fitting a straight line on log |?/^| as a function of log k, for 
the interval k = 100... 10^ we obtained 0.2554 ± 0.0006 
for /3. However for the interval k = 10^... 10^ we obtained 
0.24995 ± 0.00148, which comes even closer to the theo- 
retical value of 1/4. (II) After choosing an arbitrary initial 
point close to the focus, we moved it closer to the focus 
step-by-step by increasing k and calculating its new posi- 
tion from the old coordinates with the scaling laws for the 
fringe shrinking rate (SLG, p=l,n = — 1). In every step, 
we recorded the intensity. From the fitting on data points 
spanning the interval k = 100... 10^, we obtained that 
/3 = 0.2615 ± 0.0010, and for the interval k = 10"^... 10^ 
we got that (3 = 0.2519 ± 0.0006. The agreement with the 
theoretical exponents improves as k increases in all of the 
cases. Note that test (II) is also an indirect proof for the 
scaling laws regarding the shrinking rates. We have found 
good agreement for the scaling laws applied on the fold 
caustics as well. We note however that in BLG the wave 
pattern around the missing focus follows only the scaling 
laws regarding the shrinking rate. 

We conclude that even for the values of k we have per- 
formed our numerical calculations with, the resulting scal- 
ing exponents agree well with those predicted by the catas- 
trophe theory in the limit of /c ^ 00. 



4 Summary We have studied the scattering of a 
plane wave of electrons on a circular n-p junction in single 
and bilayer graphene with negative refractive index. We ap- 
plied the semiclassical approximation of the integral repre- 
sentation of the wave function, and tested the predictions of 
catastrophe optics on our electron-optical system. Based on 
a number of numerical calculations, we have demonstrated 
that the exact wave functions are in good agreement with 
the approximated results close to the caustics in single- 
layer graphene. Defining a quantitative error for their de- 
viations, we have found that the agreement improves in the 
short wavelength limit. We have also verified the scaling 
laws describing the asymptotic (k ^ 00) behavior of the 
intensity patterns in the vicinity of the caustics. We empha- 
size that in the case of a more complex scattering geometry, 
when exact analytical or numerical calculations are inac- 
cessible, the semiclassical approach used here still offers a 
way to obtain the electron wave intensities at the caustics. 
Therefore we think that this approach might be useful for 
design and modelling of future graphene-based electron- 
optical devices. 
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